set seed 68187

local totalSims = 1000

foreach y in economy energy immigration {
quietly {
postfile buffer effSize sigLinear sigSemipar sigFE sigRE convergedRE using Working/`y'PowerSim, replace

use Data/Pew_Models.dta, clear

 
reg abd_d unit d_congress if `y'==1

keep if e(sample)==1

bysort year: egen meanUnit = mean(unit)

local bCongress = _b[d_congress]
local bConstant = _b[_cons]
local sdError =  e(rmse)

forvalues x = 1/`totalSims' {
noisily disp `"`y' "'  `x'


forvalues effSize = 0(0.005)0.05 {
noisily disp "	effect size = "-`effSize'
gen fake = -`effSize'*unit + `bCongress'*d_congress + `bConstant' + rnormal(0, `sdError')

gen effSize = `effSize'

noisily disp "		Linear" 
reg fake unit d_congress year
gen sigLinear = ((2 * ttail(e(df_r), abs(_b[unit]/_se[unit])))<0.05)

noisily disp "		Semiparametric"
semipar fake unit d_congress, nonpar(year) nograph 
gen sigSemipar = ((2 * ttail(e(df_r), abs(_b[unit]/_se[unit])))<0.05)

noisily disp "		Fixed-effects"
reg fake unit d_congress i.year
gen sigFE = ((2 * ttail(e(df_r), abs(_b[unit]/_se[unit])))<0.05)

noisily disp "		Random-effects"
mixed fake unit d_congress meanUnit || year:,  iterate(100)
gen sigRE = ((2*normal(-abs(_b[unit]/ _se[unit])))<0.05)
gen convergedRE = e(converged) 

mean effSize sigLinear sigSemipar sigFE sigRE convergedRE

post buffer (_b[effSize]) (_b[sigLinear]) (_b[sigSemipar]) (_b[sigFE]) (_b[sigRE]) (_b[convergedRE]) 

drop fake effSize sigLinear sigSemipar sigFE sigRE convergedRE
}
}

postclose buffer

}
}

use Working/economyPowerSim, clear
gen issue = "Economy"
save Working/economyPowerSim, replace

use Working/energyPowerSim, clear
gen issue = "Energy"
save Working/energyPowerSim, replace

use Working/immigrationPowerSim, clear
gen issue = "Immigration"
save Working/immigrationPowerSim, replace

use Working/economyPowerSim, clear
append using Working/energyPowerSim
append using Working/immigrationPowerSim

collapse sigLinear sigSemipar sigFE sigRE, by(issue effSize)


set scheme s1mono

gen actualEffectSize = .
replace actualEffectSize = .01298973 if issue == "Economy"
replace actualEffectSize = .01721786 if issue == "Energy"
replace actualEffectSize = .01691328 if issue == "Immigration"

egen tag = tag(issue)

gen one = 1 if tag==1

gen point9 = .9 if tag==1

gen label = "{it:Observed}"  if tag==1

twoway ///
(line sigLinear effSize, lpattern(solid)) ///
(line sigSemipar effSize, lpattern(vshortdash)) ///
(line sigFE effSize, lpattern(longdash)) ///
(line sigRE effSize, lpattern(dash_dot)) ///
(spike one actualEffectSize, lcolor(gs9) lwidth(vthin) lpattern(dash) base(-.026)) ///
(scatter point9 actualEffectSize, mcolor(none) mlabel(label) mlabposition(3) mlabgap(0)) ///
, yscale(range(-.026 1.025)) xscale(range(-.0025 .0525)) ///
 by(issue, rows(1) note("")) subtitle(, nobox fcolor()) ///
xtitle(Absolute effect size) xlabel(0(.01).05, labsize(vsmall)) ///
ytitle("Statistical power (5% level)") ylabel(0(.2)1, labsize(vsmall) angle(horizontal))  ///
plotregion(lcolor(none) margin(0 0 0 0)) xsize(16) ysize(12) yscale(noextend) xscale(noextend) ///
legend(order(1 "Linear" 2 "Semiparametric" 3 "Year fixed-effects" 4 "HLM")  symxsize(9) rows(1) region(lcolor(none) margin(zero)))

graph export Figures/Figure_A4.pdf, as(pdf) replace
